#include "feLuSolver.h"
#include "feMvOps.h"
#include "feObserver.h"
#include "feMatrix.h"
#include "feVector.h"

feLuSolver::feLuSolver(void)
{
	A = 0;
	F = 0;
	X0 = 0;
	X1 = 0;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
void feLuSolver::BindData(feMatrix *iM, feVector *iF, feVector *iX0)
{
	A = iM;
	F = iF;

	//allocate x
	if(!X0)
	{
		X0 = new feVector(A->GetSize());
		X1 = new feVector(A->GetSize());
	}
	else if( X0->GetSize() != A->GetSize() )
	{
		X0->Free();
		X0->Allocate(A->GetSize());
		X1->Free();
		X1->Allocate(A->GetSize());
	}

	Solved = false;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
void feLuSolver::SetX0(feVector *iX0)
{
	return;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
void feLuSolver::Iterate()
{
	if(!A) return;


	//decomposition
	if(!A->LuDecomposition())
	{
		feObserverManager::GetInstance().Notify("Lu decomposition did not implemented for this matrix type");
		return;
	}

	//forward trace
	if(!A->ForwardTrace(F, X0))
	{
		feObserverManager::GetInstance().Notify("Forward trace error");
		return;
	}

	//backward trace
	if(!A->BackwardTrace(X0, X1))
	{
		feObserverManager::GetInstance().Notify("Backward trace error");
		return;
	}


	//set flag
	Solved = true;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
double feLuSolver::GetEps()
{
	if(!Solved) return double(1e10);
	return 0;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
void feLuSolver::GetX(feVector *oX)
{
	feMvOps::CopyVV(X1, oX);
}


//---------------------------------------------------------------------------------------------------------------------------------------------
feLuSolver::~feLuSolver(void)
{
	if(X0) delete X0;
	if(X1) delete X1;
}
